Printing Optimization of 3D Structure with Lard-like Texture Using a Beeswax-Based Oleogels

In this study, we investigated the optimal conditions for 3D structure printing of alternative fats that have the textural properties of lard using beeswax (BW)-based oleogel by a statistical analysis. Products printed with over 15% BW oleogel at 50% and 75% infill level (IL) showed high printing accuracy with the lowest dimensional printing deviation for the designed model. The hardness, cohesion, and adhesion of printed samples were influenced by BW concentration and infill level. For multi-response optimization, fixed target values (hardness, adhesiveness, and cohesiveness) were applied with lard printed at 75% IL. The preparation parameters obtained as a result of multiple reaction prediction were 58.9% IL and 16.0% BW, and printing with this oleogel achieved fixed target values similar to those of lard. In conclusion, our study shows that 3D printing based on the BW oleogel system produces complex internal structures that allow adjustment of the textural properties of the printed samples, and BW oleogels could potentially serve as an excellent replacement for fat.

stability and strength of the structure by changing the mass and void fraction deposited in the printed inner structure [23,24]. In this context, complex, digital food design that cannot be achieved with traditional food processing methods can become a reality [5]. However, research on the effect of infill structure on the texture of printed foods is still limited, and current research is limited to studies on 3D printed cheese [25] and mashed potatoes [5].
Therefore, in this study we sought to establish the optimal conditions for 3D food printing of alternative fats that have the textural properties of lard using alternative fats composed of BW and high oleic sunflower oil (HOSO). To this end, 3D printed objects were prepared according to the BW concentration and printing filling level of alternative fat, and the textural properties (hardness, cohesiveness, and adhesiveness) were compared with the printed lard structure.

Materials
HOSO and lard (Cypro lard oil, Daekyoung O&T, Korea) were purchased from a local grocery store. BW (foodgrade, acacia type), a yellow solid cube with a mild odor, was provided by Honest Honey (Korea).

Production of Oleogels
For the preparation of oleogels, 30 g of HOSO was weighed in a 100 ml Pyrex beaker, and an appropriate amount of BW was added according to the treatment. The mixture was heated to 90 ± 5 o C for 30 min under magnetic stirring at 300 rpm. After complete dissolution of the oleogelator, the molten oleogel was poured into a polypropylene syringe (2.2 cm inner diameter). Samples were used after being stored at 23 ± 2 o C for at least 24 h [26].

Experimental Design
A central composite face-centered (CCFC) design was used in the response surface methodology (RSM) method for the analysis of two independent variable effects (BW concentration and IL), interactions, and quadratic terms on the main properties of 3D printed objects [27]. A total of 13 experiments were performed in this work, including 5 replicates at the central point ( Table 1). The total number experiments were calculated using the following equation [28]: where N is the total number of experiments required; n is the number of factors; and n c is center runs (five replicates). Previous studies have suggested adding at least 4 center runs for a central composite design, to achieve adequate pure error freedom and reasonably sensitive lack-of-fit testing [29,30].
The sequence of experiments was completely randomized to minimize unexpected errors in responses due to external factors. BW concentration (11, 15, and 19%) and IL (50, 75, and 100%) were determined from preliminary experimental results (data not shown). Each level was selected by conducting several test experiments taking into account the maximum operating level of the printing equipment.
The dependence of each experimental response on the independent variable was modeled by applying the following second-order polynomial Eq. (2) [28,30]: where Y is the response; X i and X j are variables (i and j range from 1 to k); β 0 is the constant term; β j , β jj , and β ij are the interaction coefficients of the linear, quadratic, and second-order terms, respectively; k is the number of independent parameters (k = 2 in this study); and e i is the error [31,32].

Modeling of 3D Design
The model for the 3D structure printing of fats was designed with computer-aided design (CAD) software (SOLIDWORKS 2016) and exported as .stl files. The 3D model was designed as a rectangular prism with the following dimensions: length = 15 mm, width = 15 mm, and height = 10 mm [33]. Next, the generated .stl file was imported into the computer slicing software Simplify 3D (ver. 4.1.1). The internal fill pattern was a grid, and the internal fill density was designed with three different percentages (50, 75, and 100%) (Fig. 1). Based on preliminary experiments, the minimum percentage of the infill level was set to 50%.

3D Printing Parameters
The printing process was performed using a model 3D printer (Shinnove, Hangzhou Shiyin Technology, China) with a metal syringe connected to a temperature-controlled heating system with the temperature between 23 o C and 100 o C, and an interchangeable metal nozzle [33]. The operating parameters of the 3D printer are shown in Table 2. All printing tests were performed using printing parameters of the same value with only temperature parameters being modified to ensure the adequate printability of the sample [14]. All samples were printed at room temperature (23 o C).
Printing accuracy was calculated as the difference between the CAD design and experimental dimensions using vernier caliper (Digital Caliper, RUITOOL, China). The designed model was a square, with the same length and width (15 mm × 15 mm). The length and width were measured at three different positions, and the height (10 mm) was measured at five positions at the center and the edge. The dimensional printing deviation in each direction was determined as follows: Dimensional deviation in height (%) = ; Dimensional deviation in length and width (%) = ; where A is the height of the printed sample, B is the height of the designed model, C is the length and width of the printed sample, and D is the length and width of the designed model. Dimensions and weight values of printed samples were recorded after storage at -4 o C for 24 h.

Texture Properties of 3D Printed Products
Sample texture was determined by a texture profile analysis (TPA) of double-cycle compression test using a texture analyzer (TA.XT2i Plus, Stable Microsystems, UK) equipped with a 50 N load cell [34]. Samples in solid form and stored for 1 h at room temperature were placed on a test platform [35]. Samples were compressed every cycle to deform 60% of their original height using a 25 mm square probe at a pre-test and test speed of 1 mm/s. The hold time between compressions was set to 5 s. Next, the probe was returned to its initial position at a post-test speed of 5 mm/s. Hardness (HA), cohesiveness (CO), and adhesiveness (AD) were recorded on the test curve [33]. HA was determined as the maximum force measured during the first compression cycle in the force-time curve [34]. CO is defined as the ratio of the positive force domain (A3/A1) in the first and second compression cycles [36]. AD is defined as the area of negative force required to detach the compression plunger from the sample in the first compression cycle [37]. Each test was repeated at least five times for each type of sample.

Statistical Analysis
All results are presented as mean ± SD values. The analysis was performed in triplicate (n=3). One-way ANOVA with Duncan's test was carried out by SPSS (SPSS, USA). The graph was expressed using Sigma Plot software (version 12.5, USA). In addition, the experimental results of CCFC (central composite face-centered design) were analyzed by Minitab ® 20 software (Minitab Inc., USA). Response surface analysis modeling analyzed quadratic mathematical models including linear, squared, and interaction terms, and statistical analysis of quadratic models used analysis of variance (ANOVA). The experimental data were evaluated by various descriptive statistical analyses such as R 2 (goodness of fit of the regression model), F-value (statistical significance of the overall model), p-value, mean sum of squares, degrees of freedom (DF), and sum of squares. When the p-value was <0.05, it was considered statistically significant, and the lack-of-fit test was used to evaluate the validity of the model. The resulting data were plotted as response surface and contour plots to illustrate the relationship between the response and experimental levels of each variable used in this study.

Setting the Printing Condition Range of Variables
Prior to using an experimental statistical method for optimization, the printing condition range of variables was set. Print deviations of height, width, and length dimensions of the printed sample were calculated as the difference between the design dimensions and experimental dimensions measured using a vernier caliper to Deviation of the height dimension of printed samples ranged from -14.12 to 3.53%, which seemed to marginally affect the visual appearance ( Figs. 2A and 3). In addition, since most of the samples had negative values, a "thinner" tendency for decreases in height relative to the designed model was observed, which may be related to compressive deformation of the printed sample due to gravity [5]. BW density was observed to affect the height deviation, and generally, as the BW density increases, the deviation value decreases, indicating that printing accuracy increases. This observation is explained by the elastic rheological behavior of lard and high concentration BW oleogels. Deviation of the height dimension of lard was the smallest among the BW oleogel samples, showing the highest printing accuracy similar to that of BW-19. In the case of BW-7, -11, and -15 samples, it was observed that as the IL increases, the deviation value decreases and the printing accuracy increases, which may be due to the formation of a sturdy structure that supports the printed sample to prevent collapse as the IL increases. In contrast, the lard, BW-19 and -23 samples exhibited positive values at an IL of 100%, indicating an expanded shape of the object, which may be related to the 'die swell behavior' effect due to the robust viscoelastic properties of the above samples. This behavior reflects the fact that as viscoelastic fluid is extruded from the nozzle, the binding force on the tip wall is removed and can expand to a diameter larger than the nozzle diameter [38].
In the length and width dimensional characteristics of the printed samples (Fig. 2B), the deviations in length and width of these samples showed a tendency for 'fatter, ' somewhat positive values, which were very similar to the shape of the designed model. The highest deviations in length and width were observed in BW-7, which can be accounted for by the high spreadability due to the viscous rheological properties of this sample. In addition, as the BW concentration increased, the overall diameter variation increased, which can be explained by the "die expansion shape" of the sample with high viscoelasticity, such as height variation. Therefore, BW-15 showed the lowest deviation value of the diameter and the highest printing accuracy. In general, no significant effect of IL on the deviation of length and width was noticed. However, it was observed that in the cases of lard, BW-19, and -23, as the IL decreased, the deviation value decreased, and the printing accuracy increased.
Overall, the lard, BW-15, -19, and -23 samples were similar to the designed virtual model, indicating that food 3D printing can achieve relative accuracy with personalized designs (Figs. 2, 3).

Model Fitting of Texture Parameters
A CCFC experimental design was adopted to select the most optimal conditions for 3D structure printing with lard textural properties (HA, CO, and AD) using BW oleogel alternative fat. Before the CCFC, a pre-experiment was performed with regard to BW concentration and IL in order to set the variables (Fig. 2). Condition combinations and ranges of variables were developed to conduct the CCFC (Table 1).
In the HA and AD models, the linear term of the BW concentration (X 1 ) and IL (X 2 ) together with the interactions BW × IL were shown to be highly significant factors (p < 0.001). On the one hand, the CO model was observed as a highly significant factor only in the linear terms of BW and IL (p < 0.001). The Eqs. (3-5) models effectively account for the changes in HA, CO, and AD as a function of BW concentration and IL.  Fig. 4 visually shows the effect of BW concentrations and IL on HA, CO, and AD of the printed products. Within the studied range, HA was observed to increase rapidly with increasing BW concentration and slowly increase with increasing fill level (Fig. 4A). This result was expected because the elastic rheological properties increased with increasing BW concentration (data not shown). Additionally, an increase in IL was expected because it reflects a larger amount of oleogel being extruded to fill the internal structure of the sample. The influences of material concentrations and printing filling level on the hardness of the printed samples were similar to results reported in previous studies [5,34]. In addition, the change in hardness increased as both the BW concentration and IL increased, indicating that the hardness was influenced by the interaction of these two variables, as shown in the Eq. (3) model.

Table 3. Analysis of variance (ANOVA) for fitted quadratic polynomial models for hardness (HA), cohesiveness (CO), and adhesiveness (AD).
In the case of cohesiveness, it was observed that cohesiveness decreased as the BW increased, contrary to HA (Fig. 4B). Cohesiveness tends to decrease with increasing hardness (proportional to A1) in the ratio of the second compression cycle (A3) to the first compression cycle (A1). It can be reasoned that the harder BW oleogels initially had a more resilient structure, but once compressed they became more irreversible deformations [5].
Adhesiveness is the negative force required for the sample to separate from the compression plunger, and high negative values indicate strong adhesion [37]. A linear increase in the negative values of adhesiveness was observed with increasing BW concentration and IL (Fig. 4C), which may result because it is associated with an increased complex viscosity (rheological parameter) with increasing BW concentration [39]. The increase in infill levels was expected as an increase in oleogel extrusion volume. In addition, the adhesiveness strength obtained its maximum value when both the BW concentration and IL were increased, and as shown in the model equation, it was found that adhesion was influenced by the interaction of these two variables.

Optimization and Validation of BW-Based Oleogel Alternative Fat 3D Printing Structure
To produce optimal BW oleogels with textural properties (HA, CO, and AD) similar to lard, multiple reaction optimization was performed to determine a set of satisfactory preparation parameters that meet all the demands imposed on the reaction parameters [5]. The target values for HA, CO, and AD were fixed at 68.53, 0.47, and -146.03, respectively, which are the lard values at the 75% infill level, such that the preparation parameters and responses can take arbitrary values within the analyzed range. As a result of multiple reaction examples, the best set values were found to be BW = 15.99%, IL = 58.89%, and the according D value was 0.9945 (Fig. 5). These optimal requirements were used in texture analysis to validate the predictive model.
As shown in Table 4, the properties of printed BW oleogels obtained with the optimized parameters gave almost identical values of HA and CO from lard with no statistically significant difference. In addition, the AD was slightly higher, but was within the range of the predictive model, and there was no statistically significant difference. The BW oleogel was able to imitate the textural properties of lard, and the predictive model satisfactorily described the actual behavior.
Based on the BW oleogel composition and the printing parameter IL, we performed 3D printing optimization of alternative fat with lard textural properties (HA, CO, and AD). We found that 3D printing was able to successfully provide the designed porous structure using the BW-15 sample. The dimensional properties of lard printed with ILs of 50 and 75% and the BW-15 sample were designed with similar accuracy. The textural properties were affected by BW concentration and IL, with the BW concentration imparting the highest contribution. As the BW increased, HA and AD were found to increase in an absolute curve as opposed to CO. As IL increased, HA, CO, and AD were found to increase in a curve and then decrease slightly. In addition, the interaction of BW and IL had a noteworthy influence on all responses except CO. The BW oleogel printed with an optimization preparation variable showed no significant difference in HA, CO, and AD compared to the properties of lard printed with 75% IL. Thus, we concluded that the goal of multiple reaction optimization was substantially achieved. Additionally, the multi-response mathematical model obtained through this study enables custom printing of BW oleogels. We emphasize that the BW concentration and printing parameters can modify the textural properties of printed products, and therefore this printed BW oleogel has the potential to replace lard, which is expected to provide similar texture and sensory properties of fat for materials in need of fat replacement. Future work will carry out complex marbling-related studies with the use of protein materials based on real meat.